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Abstract 

In this work we discuss the summation of the Parquet class of diagrams within Green's 
function theory as a possible framework for ah initio nuclear structure calculations. The 
theory is presented and some numerical details are discussed, in particular the approxi- 
mations employed. We apply the Parquet method to a simple model, and compare our 
results with those from an exact solution. The main conlusion is that even at the level of 
approximation presented here, the results shows good agreement with other comparable 
ah initio approaches. 
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1. Introduction 

In recent years, considerable progress has been made in the field of nuclear structure 
calculations. Due to the increasing computational power available with better and faster 
super-computers, the regime where ah initi^ calculations are possible has increased 
considerably. For lighter nuclei, several approaches have been successful, including the 
coupled-cluster method [11-0] j variational approaches P, large-scale diagonalization 
techniques (including the no-core shell model) , the effective interactions hyperspher- 



ical method [ll| and the Green's function Monte Carlo method |12h14| . For nuclei with 
A < 4, the Faddeev/Faddeev-Yakubovsky equations flB'l provide almost exact numerical 
solutions, and such calculation are used as a benchmark for testing other approaches. 
The Green's function Monte Carlo method provides an almost exact reproduction of 
the Faddeev results flil], but the scaling of the computational effort to the number of 
particles involved makes this approach impractical for medium-sized nuclei, the current 
limit being at about A < 12. The coupled-cluster (CC) method also shows excellent 
agreement with the Faddeev/Faddeev-Yakubovsky results for '^H and *Hc [6], and this 
method scales much more favorably, currently being applied to nuclei with A < 48 [Tj]. 
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Another method with strong potential for ab initio calculations of medium-mass nu- 
clei is the self-consistent Green's functions approach 17 , . This method has a number 
of interesting features. It is possible to obtain better accuracy with a smaller numerical 
effort when compared to large-scale diagonalization approaches, and the self-consistency 
requirement provides a path to ensure the conservation of basic macroscopic quantities. 
The main advantage, however, is the close connection with experimentally available data 
sets through the elementary basic blocks of the theory, namely the many-body propa- 
gators. These contain information on several excitation processes. The single-particle 
spectral function (and spectrosco pic factor) can be experimentally extracted, for example 
in (e, e') knock-out reactions 1^ 23|, see also Refs. 21-23[. Baym and Kadanoff [24, 25| 



showed that the self-consistency requirement ensures the conservation of basic quantities 
like number of particles, energy, momentum and angular momentum. Recently, con- 
verged results for the ground state energy of "^He l 56| within the self-consistent Green's 
function framework employed by Barbieri etaZ.|17l. llSl. |53| have been presented. These 
results show that the required level of precision to meet the current standard is possible 
within implementations of Green's function based methods. 

In this paper we explore the performance of the Parquet method of re-summation of 
diagrams within the Green's functions theory to calculate the ground state energy in a 
simplified model of the nucleus. We compare our results with an exact diagonalization. 
The aim of the paper is to ascertain the possibilty of establishing the Parquet summation 
method as an alternative ab initio method for nuclear structure calculations. 

The Parquet method of summing diagrams has been known for more than 50 years, 
having first been developed by Diatlov, Sudhakov and Ter-Martirosian [26J as an aid to 
describe meson-meson scattering in particle physics. These equations have since been 
used somewhat, most notably in one- and two-dimensional electron gas calculations 27|, 



28[ . They have also been used for some critical-phenomena calculation [29 . i30j . The 



most exhaustive theoretical investigations were carried out by Jackson, Lande and Smith 
[3ll - l33 |. More recently, Yasuda has used the Parquet diagram method to construct 
approximations to the reduced density matrix of general quantum systems [sH ] . 

The effective interaction generated by this method includes a large class of diagrams, 
and requires no initial assumptions on the underlying interaction with respect to range 
and strength, as opposed to ladder type or ring type (standard random-phase approxima- 
tion) interactions. The generated interaction is symmetric, that is, the particle-particle 
part of the interaction and the particle-hole part is treated on an equal footing, thus en- 
suring that all diagrams critical to a reasonable description of the many-body system are 
included. It is therefore applicable to systems where easier approaches fail, for example 
systems where it seems that both particle-particle type (calculated by ladder/G-matrix 
approximations (36l. ISTj]) and particle- hole type (usually handled by random-phase ap- 
proximation (RPA) methods) diagrams are equally important. Systems undergoing a 
phase transition consist of one such example, and there is clear evidence that both these 
parts of the interactions play a crucial role in nuclear systems (sgI ISTj . By making a 
self-consistent calculation, the included diagram classes are all automatically summed 
to all orders. Only linked diagrams are included in the sum, which ensures that the 
method is size extensive, meaning that the total energy scales correctly with the number 
of particles .[3] . 

As mentioned above, the aim of this work is thus to demonstrate how one can sum 
the Parquet class of diagrams in a simple model and compare these results with other ab 
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initio methods such as diagonahzation methods, with a critical eye on the pros and cons 
of the method. As such, this article serves mainly as a proof-of-principle of the Parquet 
method for nuclear-physics like systems where pairing correlations and core-polarization 
play important roles. 

This work is organized as follows: in Section[5]we review some basic features of many- 
body theory and Green's function theory, before presenting our description of the Parquet 
summation method in Section [3] In Section |4] details on the numerical procedure and 
approximations are discussed. We present comparisons between our calculations and the 
exact diagonahzation solution for a simple model in Section [SJ and then we summarize 
our findings in Section [S] 



2. Green's Functions and the Interaction operator 

Parquet theory is formulated in the language of Green's function theory (or propa- 
gator formalism). This formalism is a standard framework within many-body quantum 
theory, and fuller accounts are found in most textbooks on the subject, for example in 
Refs. (38I - I40I I or [4l|. A cornparatively recent presentation is given by the text of Dickhoff 
and Van Neck from 2004 [42] ■ Here we shall be content with a short discussion of the 
one- and two-body propagators (Section 12. ip . the interaction operator (Section 12.21) . the 
self energy (Section l2.3p . and finally, the matrix inversion method of finding the one-body 
propagator (Section 12. 4p . We refer the reader to the above-mentioned references for an 
introduction to the basic concepts and for further details. 

2. 1 . Propagators 

The one-particle Green's function is defined as: 

g.0{T) = g.p{t - 1') = -z«|r{c„w4(t')}|0 

^i~z{^^\c^{t)cl{t')M) t>t' (1) 

\^{^(^\clit')c^m^) t<t\ 

Here Ca{t) and cj^{t') are the annihilation and creation operators in the Heisenberg rep- 
resentation, is the N-particle ground state and T is the time ordering operator. If 
we add a particle in state /? at a given time t' , the one-body propagator for t > t' gives 
the probability that we find the system still in its ground state if we remove a particle in 
state a at time t. Similarly, for t<t' the one-body propagator gives the probability for 
recovering the ground state when a hole is created (a particle removed) at a time t and 
then annihilated at t' . Fourier transforming to obtain the so-called Lehmann represen- 
tation, we see that the denominator is zero at energies corresponding to the excitation 
energies of the (A^ + 1) and the [N — 1) states with respect to the ground state l^*^) p^ : 

1 

9ai3{uj) = — / dTe"^''gap{T) 

^y-^^ — + y2^^ — , 

^ LV-eZ + iv ^ uj - e^. - IT] 
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where the last equation introduces the notation z"^ as abbreviation for 
(^^|ca|>I'^+^)(>I'^+^|cjj|^'|^) and so on. The energies e+ and are the energy differ- 
ences E^+^ - and - E^^^ respectively. 

The unperturbed (non-interacting, or free) one-particle propagator is given by: 

, , , ( e{a-F) 9{F-a) \ 
^ \u: — e^ + iri tj — e^ — ir]/ 

where F is the highest occupied state (at the Fermi level) in the system and e° is the 

unperturbed energy of the state \a). In this case the energy differences between the 
energy of the state with N particles and the states with A^± 1 particles is just the energy 
of the single-particle state added or removed. 

To study the effects of interactions between the single-particle states, the following 
representation of the diagonal elements of the single-particle propagator is useful 

5aa(w) = / — \ . (4) 

J„oo ^TT UJ' -UJ 

Here S{a, oj) is the spectral function, given by 

S{a,uj) = -i lim [c/c<a(w -|- ir{) - 5ac<(w - it])]. (5) 

77->-0+ 

The hole part of this is valid for energies w less than the lower Fermi energy = 
Eq — Eq~^, and is given by 

Sh{a,uj) = -lmgaa{ui) 

n 

This quantity gives the probability at a given energy lo of removing a particle (creating 
a hole) with quantum numbers a while leaving the remaining — 1 particle system at 
an energy E^^^ = Eq — lj. Similarly, the particle part Sp{a,u)) is valid for energies 
w > e+ = E^+^ - E^. It is given by 

Sp{a,uj) = --Imgaaii^) 

= E K*^+^l4l*^)|^<5(a; - (E^^ - O), 

m 

and is the probability for adding a particle with quantum numbers a to an A''-particle 
system with energy w, resulting in an -|- 1-system with energy E^^^ = Eq + w. 

For a given single-particle state, we can define the occupation number n{a) and the 
depletion number d{a) as 

n(a) = «|4c„|M'o~)= /'"do; Sh{a,uj), (8) 
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and 

= «|c„ct = / duj Sp{a,w), (9) 



respectively. It can be shown that n{a) + d{a) — 1, see for example Ref. [42|. 

In a non-interacting system, choosing the set \a) determined by the single-particle 
Hamiltonian Hq as the basis gives the hole and particle spectral functions a particularly 
simple form, being delta functions with height 1 at the energies corresponding to the 
eigenvalues of the single-particle Hamiltonian. 

In interacting systems, the spectral functions become smeared out. In principle, 
the number of poles in the propagator is infinite, giving a continuous distribution of 
probabilities for the energies of the iV ± 1 particle systems. As long as the independent- 
particle picture remains relatively correct (that is, if the interactions between the particles 
are weak), the spectral functions will have sharp peaks at clearly defined energies, which 
we then identify as single-particle states. 

The hole spectral function is relatively easy to compare to experimental data extracted 



from knock-out (e, e'p) reactions in nuclei (19|,|20j. In these experiments, an incident fast 



electron transfers a large amount of energy to a single proton inside the nucleus, sufficient 
to eject the proton, and the momentum profiles of this proton and the scattered electron 
are then measured. The most commonly extracted quantity is the so-called spectroscopic 
factor, defined as 

s^^J dp|(*r>plOP, (10) 

where Cp is a momentum state annihilation operator. In an independent-particle system 
the spectroscopic factor is either (unoccupied state) or 1 (occupied state). When the 
spectroscopic factor is less than 1, it can be thought of as measuring the amount of cor- 
relation present in the TV-particle system, being the difference between the independent- 
particle spectroscopic factor of 1 and the measured value. A word of caution is however 
necessary here. Experimentally, spectroscopic factors are defined as the ratio of the ob- 
served reaction rate with respect to the same rate calculated assuming a full occupation 
of the relevant single-particle states. They are therefore often interpreted as a measure 
of the occupancy of a specific single-particle state. However, from a strict theoretical 
point of view spectroscopic factors are not occupation numbers but a measure of what 
fraction of the full wave function can be factorized into a correlated state (often chosen 
to be a given closed-shell core) and an independent single-particle or single-hole state. 
Large deviations from the values predicted by an independent-particle model, point to a 
strongly correlated system. 

In our formalism, the spectroscopic factor is given by the height of the spectral func- 
tion at the energy of the j^*^^^) state (this follows from the orthogonality of the basis 
set \a)). 

From applying the equation of motion for a Heisenberg operator, dcn (t) / dt = —i[ca(t), H] 
to Eq. ([T]), one obtains the first step in the Martin-Schwinger hierarchy [4J], relating the 
N-|-l-particle propagator to the N-particle propagator. Thus, relating the two-particle 
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propagator to the one-particle propagator j4l|: 

= S{t - t')Sa,p + ea9al3it - t') ^^-^ 



This generates a term containing the 4-point Green's function, defined by 

= {aP\K{t^,tp-t^M)\l5). 



(12) 



Since the 4-point Green's function is antisymmetric under exchange of indices, it is 
possible to define matrix elements of K between antisymmetric two-particle states as 
shown in the last equivalence in Eq. (jl2p . provided the time arguments are exchanged at 
the same time. 

Depending on the ordering of the time arguments, the 4-point Green's function de- 
scribes the propagation of either two-particle (pp), two- hole (hh) or particle-hole (ph) 
excitations. 

The Fourier transform of K is defined as 

{al3\K{LOa,u!p,uj^,uJs)\^5) = 

and the inverse relation as 

{al3\K{tc„tp;t^,ts)\-rS) = 
1 



-\-ioo 

dcj„dLj;3dw^da;5e-'"°*°-''^^**'^+*"^*^+''^**''(a/3|i^(w„,W;3,w^,cj5)|7(5). (14) 

'too 



When the Hamiltonian H is time-independent, K does not depend on the sum of the 
time variables, meaning that it is independent of the variable t — j{ta + + + ts). 

Before proceeding, we give a short summary of our Feynman diagram rules. The 
basic building blocks are interactions, represented by a horizontal line (dotted, wavy or 
thick-lined) connected by vertical particle/hole lines. The translation of the interaction 
lines are as matrix elements, the convention for numbering the incoming and outgoing 
states always being as shown in Fig. [1] Lines are associated with the full propagator. 



<12IVI34> = / ^ 

3 4 

Figure 1: The convention for numbering the legs of an interaetion. 



implying that particle/hole lines can be translated as either particle or hole propagators, 

6 



each diagram giving rise to two Goldstone diagrams. A diagram with two parallel lines 
corresponds to two (skeleton) Goldstone diagrams, one with two particles and one with 
two holes propagating. A ring-like structure translates into a particle-hole pair. When 
combined, one of the two-particle/two- hole parallel lines might sometimes be part of a 
larger particle-hole bubble, and for aesthetic reasons no longer straight. No ambigu- 
ity should arise form this. Matrix elements of more complex operators than a simple 
interaction are represented as either rectangles (composite interactions) or circles (prop- 
agators), all conforming to the same numbering convention as the interaction matrix 
elements (Fig. [T]). In some diagrams we use arrows on the internal lines. These are in- 
tended as a graphical means of showing which states are incoming and which are outgoing 
in the matrix elements of each operator, and thus have no separate physical translation. 

2.2. The interaction operator 

Using standard many-body perturbation techniques (40|. we can obtain a diagram- 
matic expansion for the four-point propagator in Eq. ([T^ . We see that there are two 
classes of contributions, one unconnected class in which two different particle lines prop- 
agate without any interaction between each other, and a second group where the particle 
lines are connected by interaction lines. In Fig. [2] we sketch in a schematic way the two 
classes. 

The four-point interaction vertex r*"P' is called the interaction operator, defined as 
all two-line irreducible diagrams with fully renormalized propagators. To lowest order, 
p4-pt i(jgj^tical to the two-body interaction V. To make the expressions a little more 
readable, we will henceforth use roman numerals on the incoming and outgoing states, 
and Greek letters for intermediate states. We can then express the diagram for K given 
in Fig. [5] in terms of the four-point interaction vertex r^"P* as: 

{l2\K{ti,t2; t3, ^4)134) = i[gi3{ti - t3)g24(i2 - U) ~ gu{ti - t^)g2i{t2 ~ t^)] 
- dta dtp / dt^ / dt.y ^ gia{ti - ^0)52^(^2 - tp) 

•I J J J aP'jS 

X {al3\T^-P\tc,,tp;t^,ts)\jS)g^3{tj ~ t3)gs4{ts - t^). (15) 
We define the Fourier transform of the interaction operator as 

{al3\T'^~P*'{uJa,ujp,ujj,uJs\jS) = 

dto. I dtp [ dt^ [ dtse''^-'''e"^^'^e-"^-''-'e-''^'''{aP\r'^-P\t„,tf},t^,ts)\-/S). (16) 



If the bare interaction is time- independent, it does not depend on the energy, and con- 
sequently the interaction operator conserves energy and depends only on the incoming 
energies. 

The Fourier transform of the four-point Green's function can be written as 

(121/^(^1, a;^,a;3,a;4)|34) = 2TTiS{uji + uj2 - 1^3 - W4) 

X [27r(5(wi - a;3)gi3(cji)g24(i^2) - 27r(5(wi - W4)gi4('^i)g23(i^2)] 

- ^ gia{(^l)g2l3{i^2){aP\^'^~^\(^a,i^l3,UJ.y,UJs)hS)g^3{uj3)gS4{(^A). (17) 




1 2 1 2 




+ p4-pt 
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Figure 2; The four-point Green's function K, separated into a set of unconnected diagrams ii"" and a set 
of connected diagrams. The unconnected diagrams can be summarized as consisting of two unconnected 
fully renormalized propagators and their exchange contributions, as shown in the lower equation. A line 
indicates the result of a Wick contraction, depending on the time ordering this could be either a particle 
or a hole. The arrows on the lines are meant to clarify the relationship between the matrix elements and 
the pictorial description, and do not distinguish particles from holes, any of the lines could be of either 
type. 



We call the first part of the four-point propagator (the first term in Eq. (jlSp) for 
the non-interacting or free four-point propagator. It consists of a product of two one- 
particle propagators. From the above expression we see that the Fourier transform of 
the non-interacting four-point propagator is given by: 

(12|if°(wi,a;2, W3, W4)|34) = 2niS{uji + W2 — 1^3 — W4) 

X [27r(5(wi - L03)gi3{u!i)g24{uj2) - 27r(5(wi - a;4)5i4('^i)ff23('^2)] • (18) 

2.3. Self energy 

To find an expression for the one-particle propagator, we once again use standard 
many-body perturbation techniques 3^ 4^ 4^. This gives the Dyson equation, giving 
a decomposition of the propagator in terms of the irreducible self energy E (also called 
the proper self energy or the mass operator) : 



^9lji^)^{l,S;u})gsf3iu}). 



(19) 



The self energy is the one-line irreducible diagrammatic insertions to the one-particle 
propagator, as shown in the diagrammatic representation of the Dyson equation in Fig.[3l 
By iterating on this we generate the exact one-particle propagator, provided the exact 
irreducible self-energy can be found. This is unfortunately in general not possible. 



In terms of diagrams, a one-particle propagator including self-energy insertions is 
called a dressed propagator, often drawn as a double line. In the case of our Parquet 
diagrams, however, all propagators are dressed, and for the sake of simplicity, we have 
chosen to draw them as single lines nonetheless. The only exception being in Fig. [31 
where the single line represents the unperturbed propagator. 



/\ 



+ 




Figure 3: The diagrammatic representation of the Dyson equation. The single line represent the unper- 
turbed propagator, the double line represents the full (dressed) propagator. 

We can find a useful relation between the self energy and the interaction operator F^'P* 
from the equation of motion of the one-particle propagator given in Eq. ([TT|) . Inserting 
the expression in Eq. (jlSp for the 4-point propagator, we obtain: 



d 



gi2{t - t') = 5{t - t')6i2 + eiguit - t') ~ iJ2{^a\V\Pl)9Mt ^ t^)9i2{t - t') 



a/37 



X {s^\r''-p'its,^,t^,u)\^i^)g^2it^ - t') 



(20) 



Taking the Fourier transform of the above expression and performing some algebra, see 
Ref. [42] , we arrive at an expression for the one-particle propagator which is identical to 
the Dyson equation, provided we make the identification 



E(1,2;l.) = -z 



Jet 2vr ^ 



1 f duji f duj2 

2 / ^ / "27 



{la\V\(3j)gf:ss{i^i)9i,,{i^2) 
X ((5/i|F''"P*(wi, W2, wi -I- W2 - w)|2i')5^q,(i:ji + ILJ2 - w). (21) 



Here the integral in the first expression is a contour integral along the real axis to be 
closed in the upper half plane, as indicated by the C f subscript. The expression for S 
is shown diagrammatically in Fig. 2] Equations pT)) and (I^TI) together with the Dyson 
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Figure 4: The self energy E expressed by the interaction operator F'* P* . The propagators are dressed 
propagators. 



Eq. ([T9|) give the exact description of the one-particle propagator if the single particle 
propagators in the interaction operator are dressed, that is, the self energy insertions are 
included. This gives a set of non-linear equations, and any solution procedure needs to 
include some sort of self-consistency scheme, as will be further discussed in Sections 13.41 
andm 

2.4- The eigenvalue equation method 

We can write the Dyson equation, see again Eq. (jl9p . as a matrix equation using 
the notation [g] as shorthand for the matrix with gai3 with indices a, /?. Assuming the 
unperturbed propagator to be diagonal, with the inverse given as [g^]~^ = w — [e], the 
Dyson equation can be generically written as 

[5H] = [^.l-([e] + [SH]ri, (22) 

where Cq. represents the energies of the unperturbed Hamiltonian, [e] being a diagonal 
matrix with at the diagonal entries. From this we see that the poles of the propagator 
are the roots a;^ of the equation 

([e] + [E(c.A)])|A) =c.a|A). (23) 

Recalling the Lehmann representation of the propagator, see Eq. ([5]), we can identify these 
roots as the energies of the ± 1 systems. The residue matrix [S\] of the propagator at 



the pole uj\ is given by [42|: 



[Sx] - lim {oj - oj^Mu)] = I ,uJ ^)CM = sa|A)(A|. (24) 

^^"A 1 - (A|[I]'(a;A)]|A) 

The eigenstate (A| is the corresponding left eigenstate of the operator in Eq. (j23p . The 
left and right eigenstates are assumed to be normalized according to 

(A|A) - 1. (25) 

We assume that the propagator has only simple poles, the expression for the degenerate 
case is somewhat more involved. Now we can write the propagator as 

A ^ 
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The eigenvalue Eq. (j23|) is more complicated than an ordinary eigenvalue equation, as E 
depends on the energy and has to be calculated at the unknown eigenvalue. There is no 
longer only one solution, but a set of different which can be quite large (depending on 
the number of poles in S). Physically, this means that adding or removing one particle 
from the ground state no longer leaves the system in one definite state, rather there are 
several possible states, each with its own amplitude. The sum over these are still unity. 
The independent-particle model is no longer appropriate, however, calculations show that 
at least for nuclear system, much of the single-particle strength is still concentrated in a 
single state for states close to the Fermi energy. States further away, either deep down 
in the nuclear well or high up, closer to the continuum, get smeared out and cannot 
properly be called single-particle states any more. 

From the one-body propagator it is possible to find the energy of the ground state by 



using the so-called Migdal-Galitski-Koltun sum rule [44|(see also Boffi in Ref. [45|): 



2 

A<A_F ajS A<Af 



3. Parquet theory 

The formalism presented in the previous section requires calculations of several infinite 
sums, and thus we need some procedure to handle these. The Parquet method offers an 
approximation to the interaction operator T^~p^ which includes a large, infinite subset of 
the full set of diagrams. 

We observe that there are several different possibilities for reducing this four-time 
operator down to a two-time operator. Depending on the physical system in question, 
reductions to a ladder or a ring operator has been used to include either pphh (which 
include for example 2p2h excitations, see discussion below) or ph correlations respectively. 
However, as argued by Jackson and Wettig neither of these approaches meet some 
basic requirements of a many-body theory to be certain of convergence. These authors 
have further argued that a necessary (but perhaps not sufficient) condition for any many- 
body summation of diagram to converge, is that both pp and hh ladders and ph chains 
be summed to all orders. 

A Green's function based like the Faddeev random-phase approximation of Barbieri 



and co-workers, see for example 17|, |43|, couples ladder diagrams and ring diagrams to 
all orders for the self-energy. The parquet method offers a method of doing this in a 
fairly straightforward manner and includes more complicated 2p2h correlations. How- 
ever, in order to perform such calculations, approximations are necessary, in particular 
with regard to the treatment of the energy dependence of the propagators. These ap- 
proximations are discussed below. 

Thus we first discuss the principle behind the Parquet theory, namely the different 
channels in which iterative expressions for the interaction operator can be found (Section 
13. ip . Then we need some more notation, given in Section 13.21 before we are ready to 
discuss the possible two-time reductions of the four-point propagators in Section 13.31 
In Section 13.41 the Parquet equations are given in a form suitable as starting point for 
numerical implementation. 
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3.1. Channels: Equivalent ways of building the Interaction Operator 

We can obtain iterative expressions for the interaction operator defined in Section [^T^ 
by examining it order by order. To first order, it is just tfie bare interaction. The next 
order consists of two bare interactions connected by the non-interacting propagator , 
third order is found by connecting a third interaction by another , and so on. There 
are three equivalent ways of connecting the legs of the interactions, as shown in Fig. [51 
We name the different possibilities according to the numbering shown in Fig. [T] Thus, if 
we connect the lines 1 and 2 we are in the [12] channel or particle-particle channel, while 
connecting the lines 1 and 3 or the lines 1 and 4 give the [13] channel or the [14] channel, 
respectively. These two latter channels are called the particle- hole channels. The [12], 
[13] and [14] channels are the equivalents of the Mandelstam variables s,t and u from 
relativistic quantum mechanics [401 ]. 

A diagram contributing to the interaction operator F'^'P* either can or cannot be split 
into two disconnected parts, one containing the legs 1 and 2 and the other the legs 3 
and 4 by cutting two internal lines. If this splitting is impossible, the diagram is said 
to be simple in the [12] channel, if it is possible, the diagram is called non-simple. The 
particle-particle interaction V^^ is defined as the sum over all the [12]-simple diagrams. 
It is easily seen that the full interaction operator F'^'P* is obtained by iterating over V^^, 
as shown in the first line of Fig. [51 where the dash-dot line represents an [12]-simple 
interaction. The equation for the vertex translates into the well-known Bethe-Salpeter 
equation: 

(12lr4-P*(c^i,w^,L^3,L^4)l34) = (12lVi2(c^i,c^2,W3,W4)|34) 

X (a/3lX°K,w^,u;^,aj5)l7J)(7Jlr4-P*(u;^,c^5,u;3,W4)l34). (28) 

The factor ^ stems from the symmetry of the interaction with respect to the exchange 
of indices. 

Similarly, we define the particle- hole interaction V^^ as the sum over all diagrams 
which are simple in the [13] channel, that is, all diagrams that cannot be split into 
one part containing the external lines 1 and 3, and another containing the lines 2 and 
4. The particle-hole interaction V^'^ is defined as the sum over all [14]-simple diagrams 
(diagrams which cannot be split into one part containing the external lines 1 and 4, and 
another containing the lines 2 and 3). An example diagram and the splitting of different 
components is shown in Fig. [S] 

Each of these will give the full interaction operator if we iterate as shown in Fig. El 
where the dash-dot-dot (large dot) line represents an [13]-simple ([14]-simple) interaction. 
The Bethe-Salpeter equations corresponding to these diagrams are 

(12lF4-Pt(wi,W2,c^3,W4)|34) - (12lVi3(wi,c^2,w3,c^4)l34) 




X (folif"(t^5,w„,w^,L^^)l7/3)(72lr4-P*(c^^,w2,W5,^4)|<^4), (29) 
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3 4 

Figure 5: Example diagram broken down into components. The original diagram is non-simple in the 
[12] channel and can be split into two parts by cutting two internal lines in such a manner that one 
part contains legs 1 and 2 and the other legs 3 and 4. The results are both simple in the [12] channel. 
The upper part is non-simple in the [13] channel, and can be split into two parts, one containing legs 1 
and 3', the other 2 and 4'. Both these are simple in the [12] channel. The final composite diagram is 
non-simple in the [12] channel. 



and 

(12|r4-P*K,w^,t^3,^4)|34> - (12|V"(wi,c^2,W3,W4)|34) 




X (a5|ifOK,L^5,w^,L^^)|7/3)(72|r4-P*(w^,W2,W3,^5)|3<5). (30) 

Combining the information in these three equations, we see that the diagrams of F"*"?* 
fah into four classes. One class of diagrams consists of diagrams that are simple in any 
of the three channels, that is, these diagrams cannot be cut into two separate pieces by 
cutting any two lines. The lowest-order member of this class is the bare interaction, and 
the next is of fifth order in the bare interaction, shown in Fig. [71 We call this class /. In 
our calculations we will only include the first term in this series, that is, I = V. 

Then there is the class of diagrams which are non-simple in the [12] channel, generated 
by repeated iterations of the type shown in Eq. (pS)) . We call this class the L diagrams. 
The conventional ladder diagrams is a subset of this class. In terms of the V^^ interaction, 
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Figure 6: Iterative expressions for F^'P' in the [12]-channel, the [13]-channel and the [14]-channel. The 
internal arrows determine which states are incoming and which are outgoing in the matrix elements of 
each operator. 



we have that 

(12|L(cji,w^,W3,W4)|34) = 

Then there are the classes which are made from iterations of the types in equations (j29p 



1 = 




Figure 7: The diagram class I, class of diagrams simple in all three channels. The first diagram other 
than the bare interaction is of fifth order in the interaction. We have included only the first contribution 
(the bare interaction) in our calculations. 



and pop . It is fairly easy to show that each of the diagrams in the [14] channel class have 
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an exchange counterpart in the [13] channel class. Working with antisymmetrized matrix 
elements in either channel will thus include all diagrams, and only one of the channels 
need be included in the calculation. We have chosen to work with the [13] channel class, 
and call this the R class. The diagrams summed in a standard RPA calculation result in 
a subset of this class. Expressing R in terms of V^'^ gives 

(12|i?(wi,cj^,cj3,^4)|34) = 




With these definitions we see that the [12]-simple interaction V^^ can be written as the 
sum I+R and that the [13]-simple interaction V^"^ is the sum I + L. The total interaction 
operator can be written as the sum of all three diagram classes: 

(12|r4-P*(a;i,W2,c^3,W4)|34) = (12|/K, ^2, ^3, ^4)|34) 

+ (12|L(t^i,cj2,W3,W4)|34) + (12|i?(t^i,c^2,W3,W4)|34) 

= {12\{I + L + R){0Jl,LJ2,0J3:i^4)m. (33) 

Rewriting Eqs. ([5^ and ([55]) in terms of /, L and R we obtain 
(12|L(wi,cj2,t^3,c^4)|34) = 

X {al3\K° {uja., uj^, w^, ujs)\j6){{-f6\{I + L + R){oj^,ujs, 0J3,^^4)\M). (34) 

and 

(12|i?(wi,W2,W3,^4)|34) = 

/ ^/ 1? / E(w + 

X {(3j\K°{uJc.,uj0, cj^, uJs)\aS}{S2\{I + L + R){uj^,ujs,uj3,uj4)h4:). (35) 

The three equations ([551) . and ([55]) together constitute the Parquet equations. In 
addition to these, some scheme for treating the self energy consistently has to be made. 
This is discussed in Section [331 Furthermore, some simplifications with respect to the 
energy dependence is needed, as will be discussed below. However, before doing that, we 
need to introduce some additional notation. 

3.2. Angular momentum recoupling 

To reduce the basis space, it is convenient to introduce an angular-momentum cou- 
pled basis, that is, a basis where two single particle states are coupled to total angular 

momentum J. Each channel naturally give rise to its own coupling scheme. When nec- 

n n 

essary we indicate the coupling in the matrix elements as (12|y|34)j to signify coupling 
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between states 1 and 2 to total angular momentum J and M . Other quantum numbers 
such as for example the total isospin projection can easily be added. In our harmonic 
oscillator single-particle basis, the two-particle basis states are independent of Af, giving 
a huge reduction in computational complexity. The conventional coupling order in the 
particle-hole channels is to couple incoming state to outgoing state. This is the same 
notation as in Ref. [46], see this work for an extensive discussion and examples. 

The [12] coupling scheme is the standard coupling scheme. The antisymmetrized 
two-particle J[i2] -coupled state is given by 

n 

|12) = \{nijiliSity_i)(n2j2hs2tz2)JM) ^ 

1 (l_(_l))Ji+.^2-./ ^ ^ {nminm2\JM){hmi^Sims^3^m{) 

^ 77117712 miimi2 ms\ms2 

{hmi2S2ms2\i2m2) \nijimiliSit;,i) (g) \n2j2m2l2S2tz2) , (36) 

where the {limiiSimsi\jimi) are Clebsch-Gordan coefficients. 

Angular momentum algebra gives the following relations between matrix elements 
with different coupling schemes (see for example Ref. (46j [) 

(iV|34)j = Y^{-)n+^^+J+-'' J'^ I -J^i ^, |(12|\^|34),7. (37) 

and 

{12\V\M)j ^Yi-y+^'+^+'^'y ( •{ l(12|l/|34)js (38) 



where J = \/2J + 1. Similar relations hold between the [14] channel and the other two 
channels. 

We have also found it useful to employ a matrix notation in the [13] channel that 
allows us to formulate the Parquet equations as matrix equations. We define the ^[13]- 
coupled matrix element as 

(13|y|24)j = (12|y|34)j (39) 

Employing the Jj^sj-coupled matrix elements, the equations in the [13]-channel can be 
rewritten to a form close to matrix equations. The equation for R, see Eq. (I35p becomes 

(12|i?(wi,W2, ^3,^^4)134) - (13|i?(wi,W2,t^3,^4)|24) = 

X {6^\K°{uJS,Ua,UJj,UJp)\^){^\{I + L + B:){i^^,W2,UJS,UJ4)\^)- (40) 

This expression is easily transformed into a matrix equation by a suitable transforma- 
tion of the K'^ matrix. The expressions in the [12]-channel lend themselves to such a 
formulation immediately with J[x2]-coupled matrix elements. 
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3.3. Two-time propagators 

The four-time four-point Green's function can be reduced to a two-time operator by 
either requiring the "upper" and "fower" times of the diagram be pairwise equal (^3 = t^ 
and ti = t2), or by setting the "left-hand" times and the "right-hand" times of the 
diagram equal, that is, ts = ti and t^ = ^2- The first choice gives rise to the ladder 
reduction, while the second is the basis for the RPA approach, and we call it the ring (or 
chain) reduction. 

The reduction from the four-point Green's function to different two-time operators 
corresponds to the reductions of the Bethe-Salpeter equation for relativistic spinors to 
different non- relativistic equations as the Lippmann-Schwinger-equation [47| . 

3.3.1. The ladder propagator 

The ladder two-time reduction of the four-point Green's function is defined by 

{12\gPP^''{t^t')\34)= lim lim {12\K{t,t2,t' ,t4)\M) 

= hm lim -i{^^\T[c2, (t2)ci„ (t)4„ (i')4„ iUmo) 

12— rt t4 — >t 

=i[9i3{t - t')g2i{t ~ t') - gi4{t - t')g2z{t - t')] (41) 
- j dta j dtp j dty J dtsgiait - ta)g2l3it - t/s) 

X {aP\T^-P\t^,tp,t^,ts)\jS)g-,3{t^ ~t')gs4{ts ~t'). 

The two-time reductions imply the propagation of cither two particles (if t > t' , giving 
an intermediate state ^f^^^), or two holes (if t < t', the intermediate state then being 
^N-2y Pqj. |.]-^g moment we concentrate on the free part of the propagator, that is, the 
propagation of two non- interacting particles (as seen from the Eqs. and (|55|) . 

only this part is needed to calculate the Parquet equations). We Fourier transform the 
first terms of the expression in Eq.|3T]to obtain the free ladder propagator gQ^^'\ It can be 
found by inserting the inverse of the Fourier transform of the four-point Green's function 
(Eq. (|17l) ) into the expression for the Fourier transform and taking the appropriate limit. 
Since Q^^^^ is a function of one time difference only, we expect the Fourier transform to 
be a function of one energy only. As we will see below, the relevant total energy is the 
sum of the energy of state 1 and state 2, making coupling in the [12] channel the natural 
choice for this operator, namely 

(12|e7'*''(r!)|34) = / d(ti - i3)e''''*^-*^)(12|er'''(^i " ^3)|34) = 



= / d(ti-i3)e'"^*^^*^^ lim lim {l2\K^(ti,t2MM)m 



= / d{tl _ ig)g«^(*l-*3)g-*('^l+'^2)tlgj(w3+W4)t 



diJi f duj2 f dujy, f d(jJ4 
~2^ J '2tt' J '2n J '2tt' 

2TrS{uJi +a;2 - - LL>4:)i[2TrS{uJi - uj3)9i3{^i)g24{i^2) ~ 2TrS{uJi - uj4)gi4{u>i)g23{^^2)] ■ 

(42) 
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Changing variables in Eq. (|42]) to = wi + cja, w = - ^^^^^,uj' = - ta±iia, we 
obtain 

{12\gPP''''in)\34)=i [ ^[g,3{n/2+Lu)g24{n/2-io)-g,4{n/2-u;)g23{u;/2-Lu)]. (43) 
J 2Tr 

To proceed, we insert the expression in Eq. ([2]) for the one-particle propagator. To ease 
readability, we use the abbreviation for the overlaps {'^o \'^a\^n^^){'^n~^^\<^l\^o) 
and so on. The integrals are simple contour integrals. Of the four terms, two have poles 
on the same side of the imaginary axis, and we close these on the opposite half plane 
so that they do not contribute to the integral, leaving two contributions evaluated by 
using the Residue theorem. Thus we obtain the following expression for the free ladder 
propagator Qq^^'^{VL): 



{i2\gi^'^\n)\M)^i I ^ 



27r 



+ Cij — e„ + i?/ 



(E 



'24 



n 

E 



'23 



D./2 - w - e,^ 



E 



E 



It] 



E 



E 

k 

-'24 



f7/2 + w- 



r2/2 — cj — — h] 



^l2 + uj~tl 

.m4- 



E 



2?; 



E 



-E 

kl 

14 ^^23 



r2/2 + w — ej, — z?7 

■'23 \ 

cj — — z?; - 

''13 -^24 



n/2- 



^k 



IT]' 



il — en — e 



7n 



irj' 



E 



i77 
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3.3.2. The ring propagator 

We define the particle- hole ring propagator {12\QP^{t — i')|34) as the reduction 

(12|e;P''(t-i')|34)= lim lim \{l2\K{t,t' 

-{K\cUh)cm<){^o\4it3)c2imo)] 

= lim hm [-z«|r[c2„(t')ci„(t)4^(t')4W]lO 

[3 — >t 

- «|4(i)ci(i)|vI/o^)(vI//V|4(i'),2(i')|^^)]. (45) 



The incoming and outgoing states could in principle be any states, as the operator 
sequence in the propagator ensures that only state combinations of a hole-particle or 
particle-hole type give a non-zero expectation value. 

The last term in Eq. (HSl) results from the fact that the exchange part of the propa- 
gator closes the diagrams into unconnected ground state energy diagrams, involving the 
one-body density matrix elements. We do not want to include these in the interaction 
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operator, hence the definition (which in the Uterature often is called the polarization 
propagator). The proper exchanges of the interaction operator diagrams are generated 
automatically when antisymmetrized matrix elements are employed. 

If we insert Eq. (IT5|) into the expression for Qp'^ in Eq. P5|) . we obtain Qp^ in terms 
of the interaction operator r"*"?*: 

{12\gP'\t - t')|34) = ^[.gl3(^ - i').924(t - t') - gi^{t - t')g2z{t ~ t')] 

j dta j dtp J dtj J dtsgiait-ta)g2f3it-tp) 



X {al3\T^-P\t^,tf3,t^,ts)\-fS)g^3m - t')gs4{ts - t'). (46) 



Looking at the free part of the propagator for the moment, we find the expression for 
Qq'^ as a function of energy in a similar manner as for the ladder propagator. This time 
we define n = uji - W3, uj ^ uji - ^^l^, w' = cja + ^iL^. Then 

{12\gf{n)\3A) - z I ^[gi3(^ + r!/2)524(^ - n/2)]. (47) 

We then insert the expression for the single particle propagators from Eq. into the 
above expression, with the same notation for the overlaps as for the Qpp^^ calculation. 
Due to the sign of the energy variable, the two particle-hole terms survive in this case: 

J /TT ^ u) + S2/2 — e„ + 177 ^ uj + 11/2 — e^, — irj 

X 124 ^ f24 N-| 

^ uj-n/2-e++iT] ^ Lo-n/2-eJ -iT]^ 

k— m+ 
^1 o ^9 



E ^13 ^24 -^13 ^24 f^o\ 



km, 



fe + fim + ify' + 6; - £„ + IT]' 



Note the difference in the definition of i7 between this expression and the expression for 

3.4- Self- consistent Parquet equations 

We are now ready to write down the Parquet equations in a formulation suitable for 
implementations, in our case a formulation which depends only on one energy. In doing 
this, we make the approximation that the total energy fl is the same in the case of the 
ladder and ring propagators. The energy-dependence of the interaction can be seen as 
representing an average incoming energy. 

Writing the equation for R in the [13]-coupled notation given in Eq. we see that 
the Parquet equation can be rewritten in a compact matrix formulation as 

[nn)] = [iin)] + [Lin)] + [Rin)], 

[L{n)] = [{i + Rm][gf^'^(n)][{i + R + L){n)], (49) 

[R{n)] = [{i + L){n)][gp,\n)][{i + L + Rm]. 
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The equation for the self energy, given in Eq. ((2T|) . connects the self energy with the 
interaction operator: 



E(l,2;c.) = -z/ ^^(la|^|2/3)5aMwi) 




X ((5/i|r'*"P*(a;i,a;2,a;,a;i + a;2 - a;)|2i^)5yQ,(wi + tJ2 - t^)- (50) 



When we approximate the full interaction operator r^"P' in this equation by the Parquet 
interaction operator T = I + L + R given in Eq. (|49)) . we need to check whether we do 
any double-counting or not. When I consists of only the bare interaction, Jackson et 
al. [3l| have shown that only contributions having the bare interaction F as a top rung of 
a ladder term can be included. The other terms lead to double-counting either because 
they are simply equal to an already included term, or because the diagram is equal to 
an included term with some self-energy insertion, and therefore must be excluded. Thus 
the correct propagator is the Gq^'^'^ propagator and the correct coupling order is the [12] 
coupling, resulting in the following equation for the self energy: 



The single-particle propagator is expressed by the self energy via the Dyson equation, 
see Eq. (IT^ . repeated here for easy reference 



The propagators Qq^ and Qq in the Parquet Eqs. are expressed by the amplitudes 
Zap and excitation energies e in the single-particle propagator found from solving the 
Dyson equation, creating complex dependencies which have to be solved iteratively. 

Diagrammatically, all fourth order diagrams of the Parquet contributions to the self 
energy are shown in figure. [51 All diagrams to fourth order for the ladder term are shown 
in Fig.ini and for the ring term in Fig. 1101 The propagators are fully dressed propagators. 
Iteration by iteration, the diagrams are not generated order by order in the interaction, 
but rather staggered, all diagrams to fourth order being generated after three iterations. 

We can compare the self energy diagrams to fourth order shown in Fig. |S] with the 
corresponding Goldstone diagrams by closing the diagram by a hole line. Then we find 
that all Goldstone diagrams to fourth order are generated, in addition to several higher- 
order contributions as well (we remind the reader that the propagators are dressed, that 
is, all self energy insertions are included). Thus the Parquet method include ground state 
energy diagrams to the same level of precision as a coupled-cluster calculation including 
excitation operators to the fourth order, a CCSDTQ calculation j48i]. 




+ 




g^pi^) = g°^(a;) -f ^ff"^(w)S(7,(5;a;)55/3(w). 

7(5 



(52) 
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4. Approximations in the numerical implementation 



The equation set (HHl) and (|?T|) together with the Dyson Eq. (|19p gives a solution 
to the one-body propagator il solved self-consistently. To find a viable implementation, 
we have to make some further simplifications, especially in the treatment of the energy 
variables. 

We have employed a number of approximations to the complete Parquet solution. 
The first is to truncate the infinite sum of diagrams in the class / of diagrams which are 
non-simple in all channels after the first term, keeping only the bare interaction. The 
effect of missing the corrections stemming from the missing diagrams of fifth order or 
higher is small compared to the effects of some of our other approximations. 

We have also made an approximation on the energy dependence of our interaction 
r, namely that the energy parameter 57 is the same in the calculation of L and i?, even 
though it is defined differently relative to the energies of the incoming states in the two 
cases. The parameter thus represents an average incoming energy. 

The small imaginary part irj in the propagator Eqs. (j44p and (j48p is a mathematical 
necessity to ensure that the correct poles of the propagator is included when integrating 
over Lo, and the real physics occurs in the limit of 77—^0. Ideally, this should be done 
before a numerical implementation, but we have found that the poles of the propagator 
give rise to serious convergence problems, and so the it] factor must be present also in 
the numerical implementation. The physical limit is then found as an extrapolation of 
results for different values of 77. The effect of 77 is to move the poles away from the real 
axis. This reduces the severity of the poles and smoothens the energy dependence of the 
real matrix elements in exchange for larger imaginary parts. Thus the overall effect is to 
reduce the effects of the interaction between the particles, making the results closer to 
a mean field result. This need for an extrapolation introduces an additional uncertainty 
to our results. 

Our method for handling the poles and the Dyson equation is different from the 
implementation of the Green's functions approach employed by Barbieri et. al. 17 , 3, E^l • 



They solve the Dyson equation numerically, and then use a small number of solutions 
close to the Fermi level, typically two or three. The remaining solutions are accounted 
for in an average manner. This simplifies the calculation considerably. In their approach, 
the interaction F is handled in the Faddeev random phase approximation. 

The most influential approximation is our approximate method for solving the Dyson 
equation. As discussed above, in the exact solution, the standard Hartree-Fock quasi- 
particles and quasi- holes are no longer stable single-particle states. The energy of a 
single-particle state gets smeared out over a broad range of possible e nerg ies. While it 
is possible to find such multiple solution sets, see for example Ref. i^HJ, the ensuing 



complexity makes computations within our scheme far too demanding at present. We 
have therefore opted for an approximation where we keep only the solution closest to 
the first order energy e^ °'. The first order energies are determined from the energy- 
independent first order contribution to S 

/ ^Y.^a-t\Vmg,,s{u^), (53) 
J CI- ^'^ , 



7^ 



and 



e^-;-=e^^-t-I]^-°-(a/3). (54) 
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The energy of a given single-particle state is calculated from Eq. ([23]) with S calculated 
at the first order energy, that is, the energies are the eigenvalues of the equation 



(55) 



which is now a simple, linear eigenvalue problem. The number of eigenvalues is then 
limited to N, the number of orbitals in the basis. As a consequence, the hole spectral 
function in Eq. ([B]), will be given by 



The sum over k in this equation is limited to the number of orbitals, and the energies ejT 
are the orbital energies. The spectroscopic factors will be smaller than 1, as the coupling 
between states with different orbital numbers n will give hole spectral functions which 
have some probability of having a higher energy. As each energy can be identified with 
a definite orbital, the height of the spike at that energy gives the spectroscopic factor of 
that orbital. The sum rule that the occupation and depletion numbers of a single basis 
state must sum to 1 reduces to the condition that the sum of the hole amplitude and 
particle amplitude for a given energy sums to 1, giving exactly complementary hole and 
particle spectral functions. 

The integral over the energy of the single-particle propagator in Eq. ((5T|) is solved 
by the Residue theorem. In the first order term the interaction is energy-independent, 
and so no complications occur. To perform the integration in the second term, we 
make the assumption that the residues at the poles of the ladder propagator Gq^'^''' and 
the interaction T are small compared to the residue at the pole of the single-particle 
propagator. These residues will quickly be quenched by the denominator of the single- 
particle propagator as long as they are at other energies than the energies of the single- 
particle orbitals. Due to the approximate solution of the Dyson equation, the sums over 
n+ and k— in the expression for the single particle propagator reduce to restricting the 
summations over orbitals to either over or under the Fermi level. Thus our expression 
for the self energy becomes 



E(l,2;u;)=^(la|l^|2/3)^z^^ 

aP k 

+ E MV\P^)WJn^'\w + e:+)\6^,){6^\^iu; + e:+)\2^.)z, 




(56) 



k 




+ 2E E (l«l^l/37)(/37|Gr'(^ + e'a)l'5A^>('^A*|r(c. + e^-)|2^)z,^, 



aPKF 



+ 2 E E(l"l^l/'^)(/3^l^o'"'(^ + ^L°-)l'5/^)('5A*|r(c.-l-6{-)|2^)^, 



+ 2 E E(l"l^l/'^)(/37|Gr'(^ + ^L°-)l'5A*)('5A*|r(c. + e{„°-)|2z^)^, 



(57) 
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Here z^a is the amplitude of the single-particle propagator at the first-order energy el °- . 
We have chosen this solution method to incorporate the effect of the changes in the 
spectral function into the self energy while still conserving the total number of particles. 

The implication of the fixed-energy Dyson equation approximation is that while we 
loosen the restriction that the input basis states are 'good' states and allow the single- 
particles to become linear combinations of the chosen basis set, we still assume that the 
single-particle picture is valid, that is, the system can be described as a set of (quasi- 
)particles with a discrete energy spectrum. 

We present results for two different schemes for handling the energy dependence 
of the Parquet equations. In the energy-independent scheme, a fixed starting energy 
Ein is chosen and used as the input £7 in the equations for the two-time propagators. 
Typically, Ein has to have a value well removed from the poles of the propagators to 
obtain converged results. The generated interaction V{Ein) is analogous to the interaction 
from a conventional G-matrix calculation, with some important differences. The ladder 
terms includes both particle-particle and hole-hole ladders, whereas the G-matrix only 
contains the particle-particle part. The starting energy could correspond to the energy of 
the incoming particles or to the energy difference between a particle-hole pair, depending 
on context. Like a conventional G-matrix, eventual use of the generated interaction would 
be hampered by the need to extrapolate for starting energies for which the poles of the 
propagator destabilize the solution (i.e. where the incoming particle energies correspond 
to a pole in the propagator). The inclusion of hole-hole terms imply that also negative 
starting energies will give this effect. In our approach this destabilization can be handled 
by increasing the irj parameter significantly, thus generating an interaction applicable to 
all energies. 

In the energy-dependent scheme, a (real) energy grid is set up and the propagators 
and interactions are calculated at each value of the grid. The values of r(57) used in the 
calculation of the self energy (Eq. (j57l) ') are then interpolated to the energy = w-l-e^^-. 
This scheme includes the poles in the generated interaction in a more correct manner 
than the energy-independent scheme, but the added number of poles gives additional 
convergence problems. 

5. Application to a simple model 

In this work we have chosen to test the formalism on a simple model which captures 
several basic features of nuclei, such as a discrete single-particle spectrum, strong pairing 
correlations and core-polarization. Applications to given nuclei will be presented in 
forthcoming articles. The aim here is to expose the formalism. The structures seen in 
more realistic cases are here much simpler to analyze, making it possible to discern the 
origins of the observed features of the self energy and the spectral functions to a greater 
extent, including the relative effects of the pair-correlation and the particle- hole terms in 
the interaction. Thus the insights gained in this work are of relevance to the discussion 
of the more realistic cases as well. 

We describe the model as built around a three-term Hamiltonian. In Section 15.21 
we discuss the stability of the Parquet solutions with respect to the parameter -q. In 
Section lOl we discuss the self energy, and the results for the spectral function is presented 
in Section 15.41 We compare our results with an exact diagonalization in Section 15.51 
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5.1. Description of the model 

The model has N doubly-degenerate and equally spaced single-particle levels labelled 
by n = 0, . . . , N„iax and spin a — ±1. The Fermi level defines the boundary of the "closed 
core" , as shown in the first column of Fig. [TT] We define a Hamiltonian of the system 
with three contributions, a one-body part Hq and a two-body interaction V consisting 
of two terms, as follows: 

k(T kj jkl 

(58) 

The energy of the first level is set to 0, and the energy increases by a fixed amount for 
each level. We set this fixed level spacing to 1, and the coupling constants g and / give 
the relative strengths between the level spacing and the interaction. The first term of 
the interaction has a pair structure, and can only excite two particles at a time, as shown 
in the second column of Fig. 1111 The vacant positions are holes and are drawn as open 
circles. The second term in the interaction is a pair-breaking term, as it acts between 
pairs of opposite spins, creating excitations of the type shown in the third column of 
Fig. [11] This system can be solved exactly by diagonalization, enabling us to study the 
accuracy of the Parquet summation method. We pay particular attention to the effects 
of the relative strength between the level spacing and the interaction, and to the effects of 
increasing the number of single-particle orbitals and particles. As our implementation of 
the Parquet method so far only includes two-body interactions, we can gain insights into 
the infiuence of many-body correlations beyond the two-body level. In nuclei, the pairing 
component of the interaction is known to be strong, as seen by the success of the seniority 
scheme models {54^, and thus the pairing-only model is of interest in this context. We 
know that closed-core nuclei commonly have an interaction strength of ~ 20 — 30% of 
the level spacing [H^, so we will concentrate on the span [—1,1] of interaction strengths. 
We employ dimensionless variables in the discussion of this model. 

5.2. Convergence with respect to rj 

Since our method is based on an iterative procedure, we need to investigate the 
stability of the solution and see if the results converge as the number of iterations increase. 
The rj parameter in the two-time propagators (03]) and (j48p regulates the influence of the 
pole terms, determining the stability of the iterative procedure. Increasing values of 77 
give calculations in which the effect of the poles are increasingly removed, increasing the 
imaginary parts of the results. 

For the simplest case, a two- level pairing-only model {f = 0, N = 2 and p = 2), 
the self energy becomes diagonal. The unperturbed single-particle propagator structure 
as given in Eq. ([3]) is conserved. The parameter 77 can be set to in most cases for 
energy-independent calculations within this simplest model. Most starting energies give 
convergence to machine precision within 10-15 iterations. The exception to this is starting 
energies exactly at the poles in the Gq^'^'^ and Gq'^ propagators used in the calculation of L 
and R, respectively. The convergence properties for the 77 = calculations are less stable 
when the pairing constant \g\ increases. Energy-dependent calculations have roughly the 
same convergence properties as the energy-independent case. Except for some unhappy 
cases where an exact pole in F is encountered, the results for all g values converge for 
77 = 0. Setting rj > gives faster convergence. 
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Increasing the system size changes the convergence properties, as the number of poles 
increases. For the N — 10, p = 4 system energy-independent calculations converge to 
machine precision within 25 iterations for a starting energy of Ein = — 20. The energy- 
dependent calculations exhibit an even slower convergence pattern compared with the 
N = 2 case. 

Setting / 7^ changes the interaction from a purely pair-conserving to a pair-breaking 
interaction. For simplicity, we set g = when discussing the impact of pair-breaking. 
This amounts to having a pair-breaking contribution equal in size to the pair-conserving 
contribution. The pair-breaking term introduces some new, interesting features to the 
model, most notably in the fact that the self-energy is no longer diagonal. 

For the energy-independent scheme, no convergent solutions can be found in the 
N — 2, p — 2 system for values of / > when rj — 0. This is due to the closeness of 
the first-order energies of the two levels. At the value f — 1 they become equal. Setting 
T] > remedies the instability for values of / < 0.6. Above / ^ 0.6 the solutions start to 
oscillate as a function of the number of iterations for all rj. 

In the energy-dependent scheme, convergence is good for energy mesh grids with 10-30 
points. Increasing the number of points further destabilizes the solution for some values 
of /. Then some of the mesh points hits the poles in F due to the two-time propagators. 
The 77 = calculations do not converge. Increasing rj remedies this, then an almost exact 
match between all grid sizes above 10 is observed. The convergence is good for 77 > 
when / < 0.5, as in the energy-independent case. For larger / values, the instabilities 
leads to in general poor convergence properties in the energy- dependent scheme. 

Increasing the number of levels to 10, negative values of / (attractive interaction) 
become unstable in the range values oi f ^ —0.4 to / ^ —0.6. An 77 value as large as 5 
is needed before convergence is achieved. In Fig. [THwe show the difference En — -En-i 
in a log-scale plot as a function of number of iterations. The calculations are done in the 
energy-independent scheme in the N = 10, p = 2 system for 77 = 1 (upper panel) and 
77 = 5 (lower panel). The instability seen in the N = 2 system at / > is not present any 
longer, that is, with higher 77 values we see convergent results for all positive / values up 
to 1. 

We observe that adding the imaginary component gives a more irregular convergence 
pattern. The graphs shown are representative for the general pattern in all the calcula- 
tions. 

The N = 10, p = A shows similar, but slightly less convergent patterns for positive / 
values, and there is no instability around / ~ —0.5. All negative / values are convergent 
at 7/ = 1. 

The energy-dependent calculations are more unstable. For the N = 10, p = 2 system, 
convergence could not be obtained for any value of rj when / < —0.7, and for the N = 10, 
p = 4 system, the limit for obtaining convergent results is / < —0.5. 

Generally, the number of particles seems to have a larger impact on the convergence 
properties when the pair-breaking term is included. 

5.3. Self energy 

The pair model is ideal for studying the pole structure of the self energy, as it is easy 
to discern the effects of changing different parameters. The energy dependence of the 
self energy S can be investigated both for the energy independent and energy dependent 
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cases. The first order term in Eq. ((57)) is energy-independent, giving the energies at 
which the self energy is calculated at each iteration. The energy-dependent term has 
poles stemming both from poles in the Gq^'^^ propagator and from the interaction F, 
giving a rich structure in an exact calculation. The pair-conserving interaction gives a 
diagonal self energy. 

The pair-breaking part of the interaction generates off-diagonal contributions to the 
first-order self energy, that is, (1|S|2), (2|i;|l) 7^ in the = 2, p = 2 model. This is 
illustrated in Fig. [T3] showing the self energy matrix elements obtained by an energy- 
independent calculation with L = R = and / = —0.5. The numbering of the states in 
the matrix elements is given by the level numbers shown in Fig. 1111 

The QPP^^ propagator has eight poles (the Dyson equation solution gives only two 
energies in each sum in Eq. (O), but some will probably be quenched by the small re- 
moval/addition amplitudes. Inclusion of an energy-independent F will give pole positions 
which depend on the starting energy, as the poles in the two-time propagators also have 
this dependence. The position of the poles in the self energy is further adjusted by the 
self-consistency procedure, but due to our approximate solution to the Dyson equation, 
the number of poles remain unchanged. 

Reducing the absolute value of the interaction strength parameter / reduces the 
impact of the poles drastically. The main effect of an attractive force is to lower the 
average energy of the lowest-lying level and increasing the gap between the two levels. 
This lowers the ground state energy. The repulsive force reduces the gap by pushing the 
lowest level upwards in energy. 

In the upper panel of Fig.[T4]we have shown the matrix element (1|E|1) in the N = 10, 
p — A for calculations with rj = and with -q = 1. Increasing N and p, the number of poles 
increases as expected. We know that it is necessary to have > to obtain convergence 
in the larger models. As the number of poles in the larger systems quickly becomes quite 
large, the effect of is to dampen the pole structure of the self energy to resemble an 
average self energy. When rj becomes large, all structure is lost, and the solution found 
is a mean-field solution. The solution is different from the first-order solution, however, 
as the self energies are complex, and there is a larger number of non-zero off-diagonal 
matrix elements. A convergence failure for all values of 77 indicates that no such solution 
can be found within our scheme. 

The self energy matrix elements show a characteristic separation between the two 
groups of poles. The pole structure for positive energies is due to two-particle-one-hole 
excitations, while the pole structure at negative energies is due to two-hole-one-particle 
excitations. 

We have shown diagonal matrix elements (1|I]|1) and (3|I]|3) of the -q — 1 calculation 
in the lower panel of Fig. 1141 These are representative for the self energy corresponding 
to holes ((1|S|1)) and particles ((3|I]|3)). The off-diagonal element (1|S|2) and (2|i;|l) 
are representative for all the off-diagonal elements. They have values slightly below 0, 
and the remnants of pole structure is most prominent at negative energies. The average 
energies of the diagonal elements are determined by the input energy level. The structure 
at negative energies shifts to the left with increasing level number. Only the diagonal 
elements corresponding to hole states have a distinct shape stemming from two-particle- 
one-hole propagation. The pole structures of the hole diagonal elements and particle 
diagonal elements for the ry = calculation have different average features (although both 
types show the characteristic two groups of poles seen in the ry = graph of Fig. [T^ . 
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and this is reflected in the different shapes of the rj = 1 graphs. 

An interesting feature to notice in the lower panel of Fig. [14] is the fact that the 
ofT-diagonal matrix elements are unequal, due to the off-diagonal elements of the single- 
particle propagator. Thus the self energy matrix in the Dyson equation is a non- 
symmetric matrix when the number of levels and particles have increased. 

The element in an energy-independent and an energy-dependent calculation 

with 77 = is given in the upper panel of Fig. 1151 The influence of the poles of the L 
and R contributions to F makes the energy-dependent self energy more irregular and 
the poles more prominent for positive energies than in the energy-independent case. We 
observe that the poles at negative energies are shifted to the left. In the lower panel 
we have shown the same two matrix elements after the first iteration when rj = 1 in 
both calculations. The simpler structure and generally smaller amplitudes in the energy- 
independent r/ = case give a faster damping with rj than in the energy-dependent case, 
at least partly explaining the better convergence properties. The main contribution 
stemming from two-particle-one-hole propagation is shifted towards higher energies. 

5.4- Spectral functions 

The single particle propagator has non-diagonal amplitudes, resulting in a discrete 
spectral function. The coupling between states with different orbital numbers will give 
hole spectral functions which have some probability of having an higher energy. As each 
energy can be identified with a definite orbital, the height of the spike at that energy 
gives the spectroscopic factor of that orbital. When the interaction is 0, the particles are 
independent, and the spectral function will have one spike of height 1 at the energy of 
each basis state. The hole spectral function has thus 1 at the levels chosen to be holes, 
and at the levels chosen to be particles. The effect of changing the interaction strength 
/ is to reduce the height of the most dominant spikes, giving small amplitudes also at 
the energies of the other basis states, as shown in Fig. for a = 10, p = 4 system 
with / = —0.1. The state numbers are according to the level scheme in Fig. [11] Only a 
small amount of the hole strength is distributed to the particle states at this small |/| 
value. In the same figure we have also shown the hole spectral function obtained after 15 
iterations, confirming that the process of self-consistency does not influence much on this 
function when the interaction strength is small. The particle spectral function exactly 
mirrors the hole function. 

Increasing the interaction strength to f — —0.5 gives unstable results for 77 = 0. In 
the upper panel of Fig. [T7] we compare the results for the hole spectral function in the 
A^ = 10, p = 4 system at 77 = 0.5 after the first iteration and after 15 iterations. The two 
lowest-energy (hole) states starts out having a height close to 1, and there is very little 
strength on any of the higher-lying states. During the iteration procedure, the energy 
of the lowest state is considerably lowered, but the strength remains almost the same. 
The second hole state remains at almost the same energy, but the strength is reduced 
by around 30%, and the strength of the higher states increased accordingly. Especially 
the third state obtains a significant enhancement relative to the initial height, and also 
a reduction in energy. The energy shifts become smaller for the higher states, refiecting 
that the probability of excitations and the influence of the interaction is minimal at these 
energy scales. The effect of increasing 77 is illustrated in the lower panel of Fig. [TT] where 
we have shown results after 15 iterations of calculations with i] = 0.5, 1 and 7 in the 
A" = 10, p = 4 system with / = —0.5. The 77 = 7 results have a closer resemblance to 
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a calculation with weaker interaction strength (lower absolute value of |/|). The energy 
levels are more evenly spaced and more of the strength is conserved in the lowest levels, 
illustrating the general effect of 77 as a parameter which lessens the effect of the interaction 
and forces the solution closer to a mean-field pattern. 

5.5. Comparison with exact diagonalization 

We have compared the different approximation schemes discussed in Section 0] with 
the exact ground state energies obtained from diagonalization. We compare the ground 
state energies as a function of the ratio of interaction strength g to the level spacing as 
the number of levels increase and as more particles are added. 

For convergent starting energies in an energy-independent calculation, the different 
starting energies yield minimal differences. We have chosen to work mostly with a starting 
energy of -20, as this is convergent for all values of g we have studied. All convergent 
results for other starting energies have less than 3% deviation from the results shown. 

In the upper panel of Fig. [18] we show results for the simplest case, with N = 2 and 
p = 2, as a function of coupling strength. We compare the ground state energy found from 
an exact solution, a calculation where L = i? = 0, an energy-independent calculation at 
different starting energies and extrapolated energy- dependent values. In the L = R = 
calculations we have performed a self-consistent calculation of the single-particle energies 
and binding energy based on the input interaction V alone. 

We see that the agreement between our results and the exact is very good for small 
values of the coupling strength \g\, where the independent-particle properties are dom- 
inant. For stronger negative coupling relative to the level spacing, the Parquet calcu- 
lations underbind, probably due to the error introduced in our approximate solution of 
the Dyson equation. The differences between full Parquet and the results for L = i? = 
are very small, indicating that contributions to the ground state energy are provided by 
the first- and second-order self energy diagrams in the Parquet solution of this simple 
model. The results for the starting energies Ein — —10, —5, 5 and 10 demonstrate that 
in the pair-conserving model the dependence on Ein is very small. 

The energy-dependent results agree very well with the energy-independent data in the 
range of g values where the calculations converge. To improve the convergent range, j] 
has to have a non-zero value, giving the graph in the upper panel of Fig. lTSl A calculation 
based on several values of ij and linearly interpolated to 77 = performs slightly better at 
large g values. The tendency to underbind at more negative g values persists, showing 
that the energy dependence alone is probably not enough to recover all of the missing 
binding. 

When we increase the number of levels to 10 and the number of particles to six, we 
obtain the graphs in the lower panel of Fig. [THl We see increasing differences between 
the exact values and the Parquet calculated values at stronger negative coupling. The 
extrapolated values are still better than the energy-independent values which in turn are 
slightly better than the L = R = results. The dependence on Ein is still negligible, 
which is why only the Ein = — 20 graph is shown. 

To study the effects of increasing particle number, we show in Fig. [TO] the absolute 
difference 

I ^Parquet ^Exact | 

^ ~ TgExactl 
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between our energy-dependent results and the results from the exact diagonalization for 
the N = 10 model with two, four and six particles. The odd data point at / = — 1 in the 
p — 4 graph has a value close to 7. This rapid increase is mainly due to the fact that the 
exact value at this point is close to 0, amplifying the difference at this point. We see that 
in general, the relative errors of our results reduce as the number of particles increase. 

We have also investigated the changes in correlation energy as the number of particles 
change. We show in Fig.[20]a comparison between the energy-dependent results and exact 
diagonalization results. The two latter graphs are adjusted such that the ground state 
energies at g = are equal. Thus the graphs for p = 4 and p — 6 show the correlation 
energy Ec in the system due to the increased number of particles. As the number of 
particles increase, the correlation energy changes more rapidly with increasing interaction 
strength, as could be expected. 

In the lower panel we show the relative difference in correlation energy 

I ^Parquet _ ^Exact | 
= |£:Exact| 

for the p — A and p = 6 cases relative to the p — 2. The errors have two distinct 
sources, namely the systematic errors introduced by the approximations employed in the 
solution, and the errors stemming from missing many-body correlations. If the error 
I ^Parquet _ ^Exact | g^-aigg |;]-^g gxact Correlation energy with respect to the number of 
particles, the relative difference would be independent of particle number for four and 
six particles. In the lower panel of Fig. we have plotted the relative difference. The 
graphs for four and six particles coincide in the range —0.7 < / < 0. The / > values 
show that the Parquet solution is very close to the exact, and the different behaviour 
might originate from the uncertainties in the extrapolation to r/— 5-0. 

In the pair-breaking model, calculations for the simple N — 2, p — 2 system have 
a close correspondence to the exact calculation for / < values, as shown in Fig. 
The differences between different starting energies are rather small in this model as well. 
The upper panel of the figure shows graphs for Ein = — 5 and an extrapolation based on 
extrapolation from calculations with 77 = 1,5 and 10. For / > the results destabilize 
around / = 0.3, as expected from the discussion on the convergence. The extrapolated 
results are not able to do any better than the 77 = results, and overbinds at / < in 
the same manner as the L = R = results. Energy-dependent calculations, as shown in 
the lower panel of Fig. [2Tl show in general the same behaviour as the energy-independent 
calculations. Results extrapolated from calculations with ry = 1,2,3,5 and 10 manage 
to agree nicely with the exact solution up to / = 0.5 before becoming too unstable. As 
stated in Section [521 all calculations become unstable at / = 1. It is possible to find 
convergent solutions for values of / > 1, but they are very far from the exact value. The 
exact ground state energy decreases, while the Parquet solutions increase. The failure 
is due to the reversal of the lowest-lying levels, violating our assumption that the input 
basis should not differ too radically from the expected solution. 

The results for = 10, p = 2 and N = 10, p = 4 arc shown in Fig. [511 We compare 
L = R = 0, energy-independent and energy-dependent calculation for calculations with 
the exact solution for both the N — 10, p — 2 model (upper panel) and A^ = 10, p = 4 
(lower panel). 

In the A^ = 10, p = 2 model we see that the instability around / = —0.5 seems to 
indicate a breakdown, as both the L = R = and the energy-independent calculations 
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give far too negative ground state energies when / < —0.5. In the iV = 10, p = 4 model 
the energy-independent calculations and the energy-dependent calculations are rather 
close, and slightly closer to the exact solution than the L = R — calculations. 

All values shown in Fig. [52] are converged, but the unconverged results at smaller 77 
values show oscillations between a state with larger and a state with a much lower ground 
state energy than the exact value. The iV = 10, p = 2 results shown are consistent with 
the lowest energy, while the N — 10, p = 4 results is more in line with the highest energy 
in these (large) fluctuations. This might indicate that the extrapolated results shown 
does not represent the ground state energy at 77 = properly in either system. The 
discrepancies become significant at values / < —0.5, consistent with the indications of a 
breakdown of our solutions around this interaction strength. 

It must be remembered that all values are extrapolated and carry an uncertainty 
which roughly estimated is larger than the observed differences between the three ap- 
proaches. Thus it is more prudent to look at the general trends rather than analyze the 
finer points in any great detail. 

6. Conclusion 

We have applied our implementation of the Parquet method to a simplified model of 
the nucleus, a model with a constant spacing between single-particle levels and particles 
interacting via a constant interaction. We have looked at two schemes of the interaction, 
one in which only pair-excitations are allowed, and one with a pair-breaking term the 
same size as the pair-conserving term. In both cases we have investigated the stability 
with respect to rj when varying different parameters, and compared different energy 
approximation schemes to the exact solution found by matrix diagonalization. 

We have found that the Parquet method performs generally well with the pair- 
conserving interaction. In this model the self energy is diagonal. The stability of our 
solutions depend on the parameter 77, which regulates the influence of the pole struc- 
tures in the propagators. The convergence properties with respect to rj are normally 
quite good, with some notable exceptions when a pole in the generated interaction F is 
encountered. At small system sizes, the number of poles is sufficiently small so that no 
77 is needed, but as the number of levels increases, a finite, but small 77 is needed. The 
number of particles have little impact on the convergence properties. 

The energy-dependent scheme was found to perform slightly better than the energy- 
independent scheme when compared with the exact results. Over a range of interaction 
strengths between -0.5 and 0.5 of the level spacing both schemes show good agreement 
with the exact solution, also as the number of levels and particles in the model is in- 
creased. At larger interaction strengths, our solutions underbind the systems. In the pair- 
conserving model, the effect of increasing the number of available levels is rather large. 
It is difficult to ascertain the relative importance of the systematic errors stemming from 
the limitations imposed by our approximations when implementing the Parquet method 
versus the effect of missing many-body correlations. The study of the correlation en- 
ergy with increasing particle number indicate that our solution method scale well with 
increasing particle number. 

Introducing a pair-breaking force destabilizes the solutions in the region of strong 
repulsion (/ > 0.5) in the smallest system with two levels and two particles, with a com- 
plete breakdown of the method when the first-order energies of the lowest levels become 
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equal. Increasing the number of levels in the model gives larger differences between the 
exact solution and the Parquet solutions as the absolute value of the interaction strength 
increase. Increasing the number of particles gives increasingly larger discrepancies. 

As discussed, the energy-dependent calculation method is closer to the exact solution 
in the pair-conserving systems, but this is not the case in the pair-breaking systems. To 
obtain convergence in these latter systems, rj has to be increased to such large values 
that the solution becomes a mean-field-type solution with an almost energy-independent 
self energy. The differences between results for different 77 values are small. Most of the 
correlations have been lost, and the energy- dependent results are rather far away from 
the exact solution. 

The energy-independent calculation method gives results closer to the exact solution 
than the energy- dependent scheme. Within this simple model, there is almost no depen- 
dence on the starting energy as long as this is chosen at values where the generated 
interaction F does not have any poles. 

The general conclusion is that for the larger pair-breaking systems, the Parquet 
method as implemented is only in agreement with the exact solution for —0.3 ^ / ^ 0.2. 

The main conclusion of our study of the Parquet method applied to our simple model 
is that its best area of application is when the interaction strength is rather small com- 
pared to the level spacing. The method seems to fail to find solutions close to the exact 
solution as the number of levels and particles increase. The energy-dependent scheme is 
closer to the exact solution for small systems where the number of poles are limited, but 
as the number of poles in both self energy and generated interaction F increase when the 
system size increase, this solution scheme needs large values of 77 to converge. Then too 
many details of the pole structure are lost, and the solutions found by this scheme are 
further removed from the exact solution than the energy-independent solution. 

In a larger perspective, the Parquet method can be implemented and give reasonable 
results which seems on the right track to becoming useful as an ab initio nuclear structure 
calculation method. The simple model discussed above has many factors in common with 
realistic systems. Many of the general conclusions and insights gained can be transferred 
to an application to real nuclei. It is easy to generate graphs for the self energy as a 
function of energy, and these can be used to determine the degree to which the solution 
incorporates correlations beyond a mean-field type solution. The spectral functions are 
easy to extract, as are spectroscopic factors. 

The Parquet method certainly has the possibility for attaining a high level of accu- 
racy. The angular-momentum coupling schemes allow for huge reductions in the model 
space, giving the possibility to perform much larger calculations. All Goldstone ground 
state energy diagrams to fourth order are included, a feat which in the coupled-cluster 
approach require inclusion of excitation operators to the fourth order (implying a CCS- 
DTQ calculation) [4^. The results on the ground state energy of ^He [56] show that 
the required level of precision to meet the current standard is possible within implemen- 
tations of Green's function based methods. Results for nuclei like "'He and ^^O will be 
presented in a forthcoming work (st} . 
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Figure 9: The ladder diagrams generated by the Parquet method. All contributions to fourth order are 
explicitly drawn. The propagators are dressed propagators. 
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Figure 11: Sketch of the N = 10, p = 4 model in the ground state (1), an example pair excitation, the 
only type to occur when g 7^ 0, / = (2), and an example pair-breaJcing excitation which can occur 
when / ^ (3). 
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Figure 12: The difference En — E^—i between successive iterations as a function of the number of 
iterations for energy-independent calculations with Bin = —20 for several / values, with rj = \ (upper 
panel) and t; = 5 (lower panel) in the AT = 10, p = 2 model. 
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Figure 13: Energy-independent calculation of E as a function of energy calculated for / = —0.5 in the 
N = 2, p = 2 system. The legends shows the indices of S, that is,(l|Ejl) and so forth, according to the 
levels in Fig. 1111 
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Figure 14: The upper panel shows the (1|E|1) matrix element in the N = 10, p = 4 system at / = —0.5 
for rj = and for rj = 1 after the first iteration. The lower panel shows the diagonal matrix elements 
(1|E|1} (labelled 1) and (3|S|3) (labelled 3), and the off-diagonal elements (1|S|2> (labelled 12) and 
(2|E|1) (labelled 21). The oflF-diagonal elements are not equal. 
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Figure 15; The difference between tiie matrix element (1|S|1) of the self energy of an energy-independent 
and an energy-dependent calculation after 1 iteration in the TV = 10, p = 4 system with = in the 
upper panel and with r; = 1 in the lower panel. 
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Figure 16: Spectral function in the N = 10, p = A system with / = —0.1, after 1 and after 15 iterations. 
The numbers above each bar is the level associated with that energy (due to the approximation of the 
Dyson equation, such a correspondence can be made), according to the numbers of Fig. 1111 The duplicate 
numbers 3i and 3 and so forth indicate the spectral function after 1 iteration and after 15 iterations, 
respectively. Only a small amount of the hole strength is distributed to the particle states at this small 
I /I value. 
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Figure 17: Spectral function in the N = 10, p = 4 system with / = —0.5. In the upper panel the 
differences between 1 and after 15 iterations is shown. The lower panel shows several rj values after 15 
iterations. The state numbers are according to the level scheme in Fig. 1111 The duplicate numbers li 
and 1 and so forth indicate the spectral function after 1 iteration and after 15 iterations, respectively. 
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Figure 18: Comparison between the exact solution and the different Parquet solutions in the N = 2, p = 2 
model (upper panel) and in the A'^ = 10, p = 2 model (lower panel). The differences between the exact 
solution and the various Parquet solutions are much larger in the N = 10 case. The energy-dependent 
results are slightly closer to the exact solution than the energy-independent results. 
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Figure 19: The relative diflference in energy e = |EParquet _gExact|/|£;Exact| for p = 2, p = 4 and p = 6. 
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Figure 20: Comparison between the iV = 10, p = 2 and adjusted p = 4, p = 6 exa<;t and Parquet 
energy-dependent results are shown in the upper panel. The lower panel shows the relative difference in 
correlation energy sc = | eP^'''''"''* - E^'">^^\/\E^''^''^\ for p = 4 and p = 6. 
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Figure 21; Ground state energy of calculations in the pair-breaking model compared with exact solution 
for the A' = 2, p = 2 system. The results in the upper panel are for energy-independent calculations, 
while the lower panel shows energy-dependent results. For negative values of /, the agreement between 
our calculations and the exact solution is very good. Our calculations become increasingly unstable for 
positive values of /. ^ 
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Figure 22: Ground state energy of calculations compared with exact solution for the TV = 10, p = 2 model 
(upper panel) and = 10, p = 4 (lower panel). The results shown are convergent, but unconverged 
results for lower rj values indicate that the extrapolations probably do not reflect the correct 77— ^-O limit. 
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